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Abstract 


Recently, a parallel pathway model to describe ankle dynamics was proposed. 
This model provides a relationship between ankle angle and net ankle torque as 
the sum of a linear and nonlinear contribution. A technique to identify parame- 
ters of this model in discrete-time has been developed. However, these param- 
eters are a nonlinear combination of the continuous-time physiology, making 
insight into the underlying physiology impossible. The stable and accurate esti- 
mation of continuous-time parameters is critical for accurate disease modeling, 
clinical diagnosis, robotic control strategies, development of optimal exercise 
protocols for longterm space exploration, sports medicine, etc. 

This paper explores the development of a system identification technique 
to estimate the continuous-time parameters of ankle dynamics. The effective- 
ness of this approach is assessed via simulation of a continuous-time model of 
ankle dynamics with typical parameters found in clinical studies. The results 
show that although this technique improves estimates, it does not provide ro- 
bust estimates of continuous-time parameters of ankle dynamics. Due to this 
we conclude that alternative modeling strategies and more advanced estimation 
techniques be considered for future work. 
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Nomenclature 


Acronyms 

IRF- Impulse Response Function 

NARMAX- Nonlinear AutoRegressive, Moving Average exogenous 

ELS - Extended Least-Squares 

PEI - Prediction Error Identification 

SLS - Standardized Least-Squares 

PRBS - Pseudo-Random Binary Sequence 

SNR - Signal-to-Noise Ratio 

NF - Noise-Free 

CT - Continuous-Time 

FIR - Finite Impulse Response 

HR - Infinite Impulse Response 

LTI - Linear Time Invariant 

Symbols 

, ! X l - nonlinear mapping 
u - exogenous input 
Z - output 

e - innovation, or uncontrolled input 

n - sample index point 

1 - maximum polynomial order 

n z - maximum output lag 

n u - maximum input lag 

n e - maximum error (innovation) lag 

I- inertia 

B - viscosity 

K - elasticity 

A - reflex time delay 

C - damping parameter 

c 0 - natural frequency 

g - gain 

Co - zeroth-order coefficient of polynomial nonlinearity 
ci - first-order term coefficient of polynomial nonlinearity 
C 2 - second-order term coefficient of polynomial nonlinearity 
s - Laplace variable 
U ( s ) - ankle angle 

V ( 5 ) - ankle velocity 

V (s)e A ' - delayed ankle velocity 
X ( s ) - output of half-wave rectifier 
Yl(s ) - noise-free linear path output 
Fv/.(V) - noise-free nonlinear path output 

Y ( s ) - noise-free net torque 
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E (s) - measurement noise 

0 O — ^ 11 - discrete-time coefficients of ankle dynamics model 
T - discrete-time delay 

V ( n ) - terms of discrete-time ankle model associated with nonlinear polynomial 
coefficient c\ 

X{n) - terms of discrete-time ankle model associated with nonlinear polynomial 
coefficient ci 
T - sampling rate 
£ - prediction errors 
M - set of real numbers 
N - data length 

Z - vector of measured outputs 
Z - vector of predicted outputs 

0 - extended least-squares estimate of the system parameters 

p - number of system parameters 

V F - partitioned regressor matrix 

y ¥ zu - regressor matrix that is a function of z and u only 

- regressor matrix that represents all the cross products involving e 

- regressor matrix that is a polynomial function of the prediction errors 
'P - centered and standardized variate of l F 

Z - centered and standardized variate of Z 
Zip - diagonal matrix of standard deviations 
Zip*. - standard deviation of the kth column of V P 

- matrix who’s kth column has all entries equal to the mean of column k of 

0 - centered and standardized system parameters 



Hz - Hertz 
ms - millisecond 
dB - decibel 
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1 Introduction 


In a zero or low-g environment it is critical to keep astronauts healthy and 
functional. To maintain bone and muscle mass, astronauts need to exercise ev- 
eryday [1], Under terrestrial 1-g conditions gravity works against muscles and 
bones, which requires the neuromuscular system to maintain enough muscle 
and bone mass to support body weight. In an altered gravity-inertial environ- 
ment the forces of gravity are significantly reduced. As a result, astronauts 
lose muscle mass and bone density since it is not required to support body 
weight [2], 

Under current technological limitations (e.g. lack of artificial gravity through 
inertial forces), rigorous physical activity is the only successful way to compen- 
sate for the lack of gravity [2], Nevertheless, even with rigorous exercise, as- 
tronauts have typically lost 0.4-1% of their bone density per month in space [3]. 
Upon return to Earth, with a thorough physical therapy regime astronauts pro- 
gressively recover muscle tissue and much of the bone mass lost during a semi- 
long duration (i.e. 4-6 months) in low Earth orbit. However, throughout the 
mission span it is important that astronauts are strong enough to perform stren- 
uous activities in space, such as spacewalks and emergency procedures during 
landing. Due to these possible vital mission activities, a regular exercise routine 
in orbit prepares astronaut for such situations and accelerates a reconditioning 
period to recover muscle and bone loss. 

A parallel pathway model of ankle dynamics has been proposed that mea- 
sures changes to intrinsic, reflex and muscle properties [4], These characteris- 
tics change, for example, due to a lack of effective exercise. Therefore, astro- 
nauts health and exercise effectiveness can be monitored in orbit by evaluating 
deviations of these parameters from optimal pre-orbit measures. In this work, 
we propose to use this parallel pathway model as a foundation to assess lower 
limb health during orbit to maintain an optimal exercise program. 

Traditional approaches to nonlinear system identification of human ankle 
dynamics have relied on quasi-linear methods, e.g. impulse response function 
(IRF) method [5]. These methods provide convenient, robust means of charac- 
terizing the dynamics of nonlinear systems without requiring a priori assump- 
tions regarding the system structure. However, nonparametric techniques may 
require many parameters to describe even simple systems and can be difficult 
to relate to the parameters of the underlying physiological system. 

The NARMAX (Nonlinear AutoRegressive, Moving Average exogenous) 
model structure has been shown to be well suited to modeling the input-output 
behavior of ankle dynamics [6]. The unknown model parameters can be esti- 
mated in discrete-time using the extended least-squares (ELS) algorithm [7- 
9]. However, discrete-time parameters are a nonlinear combination of the 
continuous-time physiology, making insight into the underlying mechanisms 
difficult. The primary objective of biological system identification is to provide 
parameters that give insight and relevance to the underlying biology, which are 
represented in continuous-time. 

The problem of continuous-time system identification from sampled input- 
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output data can be divided into two broad approaches: (i) indirect methods, 
where a discrete-time model is estimated from sampled data; then an equiv- 
alent continuous-time model is calculated and (ii) direct methods, where a 
continuous-time model is obtained directly without going through the inter- 
mediate step of first determining a discrete-time model; based on concepts 
of approximate numerical integration to recreate time -derivatives needed in 
continuous-time formulations [10]. 

The central problem for indirect methods is that of estimating stable and 
robust continuous-time parameters from discrete-time estimates. For both lin- 
ear and nonlinear systems, the inverse mapping problem often results in large 
bias or sign flips. To address these issues several indirect techniques have 
been developed to modify the estimators, keeping the computational complex- 
ity low [11-13]. In addition, a method for estimating the natural logarithm of 
a square matrix to map discrete-time parameters to continuous-time has been 
proposed [14]. Most other procedures are variants of these two approaches. 

Direct methods avoid the inverse mapping problem by identifying a continuous- 
time model directly from sampled data. This raises several technical issues. 
Unlike a difference equation model, a differential equation model contains 
time-derivative terms that may be required but not available from measure- 
ments. Many techniques have been devised to deal with the need to reconstruct 
these time-derivatives (see e.g. [10, 15-17]). However, these approaches are 
challenged by a lack of robust implementation for numerical computation of 
derivatives, selection of filter cutoffs to suppress high frequency noise, filter 
order, etc. These problems are more complex when the system under study is 
nonlinear, which is inherently the case for biological systems. 

We propose to investigate the problem of continuous-time system identifi- 
cation by developing methodology using an indirect approach. One such ap- 
proach relies on matrix preconditioning techniques, such as standardized least- 
squares, to improve the spectral properties of the regressor matrix [18, 19]. The 
transformed matrix will have a smaller spectral condition number, and eigen- 
values clustered around one. Often when the clustered spectrum is away from 
zero it results in rapid and robust convergence, especially when the precondi- 
tioned matrix is close to normal. We deem this will provide more robust solu- 
tions and greater consistency for nonlinear biological systems by circumvent- 
ing issues of implementing numerical derivatives and filter selection required 
by direct techniques, which are more challenging in a nonlinear framework. 
Here, we focus on one critical issue, namely, that of mapping the underlying 
continuous-time system to discrete-time for estimation, then inverse mapping 
back to continuous-time to provide physiological relevance and insight. 

The organization of this paper is as follows. The NARMAX model struc- 
ture is described in section 2. Section 3 discusses a continuous-time representa- 
tion of a parallel pathway model describing ankle dynamics and its NARMAX 
representation. In Section 4 we summarize the parameter estimation method 
implemented in this study, namely, extended least-squares. Section 5 shows 
how to standardized the identification data to formulate a standardized least- 
squares approach, which we investigate for its utility as an indirect technique 
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to estimate continuous-time parameters of ankle dynamics. Section 6 illustrates 
the results of our proposed approach on a simulated ankle model. Section 7 pro- 
vides a discussion of our findings and suggestions for future direction. Section 
8 summarizes the conclusions of our study. 


2 NARMAX Model Description 

The input-output relationship of many nonlinear dynamic systems can be writ- 
ten as the nonlinear difference equation 

z(n) = & l [z(n- l),--- ,z(n-n z ),u(n),--- ,u(n-n u ), (1) 

e(n — 1), • ■ • ,e(n — n e )] +e{n) 

where & 1 is a nonlinear mapping, u is the exogenous input, z is the output, 
and e is the innovation, or uncontrolled input. This structure describes both the 
stochastic and deterministic components of a system. This nonlinear mapping 
may include a variety of nonlinear terms, such as terms raised to an integer 
power (e.g. u 2 (n — 3)), products of present and past inputs (e.g. u(n)u(n — 1)), 
past outputs (e.g .z{n— l)z(n — 2)), or cross-terms (e.g. u 2 (n — l)z(ra — 2)). This 
system description encompasses most forms of nonlinear difference equations 
that are linear-in-the-parameters. 


3 Parallel Pathway Model of Ankle Dynamics 

The Neuromuscular Control Laboratory (Department of Biomedical Engineer- 
ing, McGill University, Montreal Canada), has developed a parallel pathway 
model (Fig. 1) to describe ankle dynamics [4], This model provides a rela- 
tionship between ankle angle (rad) and net ankle torque (Nm) as the sum of a 
linear and nonlinear contribution. The upper, linear pathway models intrinsic 
stiffness as a second-order system with parameters corresponding to inertia (/), 
viscosity (B) and elasticity (K). The lower, nonlinear pathway models reflex 
stiffness as a cascade of a derivative, a reflex time delay, a static nonlinearity 
(i.e. half-wave rectifier) and a low-pass system representing muscle activation. 
The parameters associated with the low-pass system are damping parameter 
(£), natural frequency (co) and gain (g). This model is ideally suited to study 
the effects of exercise (in space) because it quantifies both the intrinsic and re- 
flex stiffness components. The parameters of these two paths will deviate from 
normal a priori terrestrial measurements due to a suboptimal exercise regime. 
A method to monitor changes of these parameters could be used to develop 
individualized exercise protocols for astronauts. 

Kukreja et al. [6] showed that a second-order static polynomial (cq +c\x(n) + 
C 2 X 2 (n)) provided a good approximation to the half-wave rectifier which gave 


6 



Figure 1 . Continuous-time morphological model describing ankle dynamics. 

Upper (linear) path: intrinsic stiffness. Lower (nonlinear) path: reflex stiffness. 

a NARMAX representation for this model of ankle dynamics as 

z(n) = 0 o + 9 iz(n — 1 ) + 02 z(n — 2 ) + 0322(72) + O4 u(n — 1 ) + 0522(72 — 2 ) + 6(,u(n — 302 ) 
+ 0722(72 — 4 ) + 0 g [22(77 — t) + 22(22 — T — 1 ) — 72(77 — T — 2 ) — 22(77 — 1 — 3 )] 

+ 0 g [m 2 (n — t) + 3 u 2 (n — x — 1 ) + 3 m 2 (« — T — 2 ) + u 2 (77 — T — 3 ) 

— 222(72 — T)n(/2 — T — 1 ) — 4 n(n — T — 1 )u(n — T — 2 ) — 2 u(n — T — 2 )u(n — T — 3 )] 
+ 01 O <?(72 — 1 ) + 0n<?(?2 — 2 ) + 22(72) 

= 0 o + 0 iy(?i — 1 ) + 02 y(« — 2 ) + 03 22(72) + 6411(11 — 1 ) + 6511(11 — 2 ) + 6^u(n — 3 ) 
+ 0722(72 — 4 ) + 6$v(n) + 6 9 %(n) + 6[ 0 e(n - 1 ) + 0 ne (/2 - 2 ) + 22(72) 

where T is the discrete-time delay, 

v(n) = u(n — t) + 22(72 — T — 1 ) — 77(72 — T — 2 ) — 22(72 — T — 3 ) ( 3 ) 

and 

^(72) = 22 2 (72 — t) + 372 2 (?2 — t — 1 ) + 3 u 2 (n — T — 2 ) + u 2 (n — T — 3 ) ( 4 ) 

— 222(72 — t)u(ii — T — 1 ) — 422(72 — T — 1)22(72 — T — 2 ) 

— 222(72 — T — 2)22(72 — T — 3 ). 

Table 1 gives the relationships of these discrete-time NARMAX parameters 
(Eqn. 2 ) to the underlying continuous-time coefficients (the “extra” coeffi- 
cient T denotes the sampling rate). An estimate of unknown system parameters 
can be obtained using standard prediction error identification (PEI) techniques, 
such as extended least-squares. 
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NARMAX 

Relationship to 

Coefficient 

Continuous-time Coefficient 

0o 

4c 0 gftPr 2 
4+w 2 T 2 +4^coT 

0i 

-8+2 co 2 T 2 
4+a 2 T 2 +4t^(oT 

02 

- 4£(qT+4+co 2 T 2 
4+ffl 2 r 2 +4for 

03 

j2 + Y +K 

04 

t— 21 B\ K -8+2fl) 2 r 2 \( I 1 B 1 

V T 2 T> U 4+co 2 T 2 +4t^coT A 7-2 + t 

05 

f 1 \ II -8+2 io 2 T 2 \(—2I B\\ 11 -4^mT+4+a) 2 T 2 ^, / 1 B 1 t'W 

^T 1 ’ U 4+co 2 T 1 +4£g>T "T 2 " T>) U 4+a i T 1 +4£a>T 2 7 +K >> 

06 

t — 8+2o) 2 7’ 2 w 7 \ f( -4£(i>r+4+ffi) 2 7’ 2 w-2/ B\\ 

v 4+(D 2 T 2 +4^mT >yT 2 > H 4+a 2 T 2 +4^coT A T 2 T>> 

07 

K —4£a)T+4+a 2 T 2 w I \\ 
A 4+m 2 T 2 +4^mT )\T 2 >) 

08 

geo 2 T 2 ci 

(■ 4+io 2 T 2 +4£coT)T 

09 

garT 2 C 2 

{4+m 2 T 2 +4^a>T)T 2 


Table 1. Theoretical relationship of NARMAX model parameters to 
continuous-time system coefficients. 

4 Parameter Estimation 

Many parameter estimation techniques are based on least-squares theory. How- 
ever, ordinary least-squares algorithms for linear systems cannot be applied 
here because they assume that the noise terms in the model are independent 
and the regressor matrix is deterministic [20]. Both of these conditions are vi- 
olated in Eqn. 2, when one considers that only a noisy measure of the output is 
available. 

Clearly, the regressors in Eqn. 2 contain correlated error terms and noisy 
data. To obtain unbiased parameters other estimation techniques based on least- 
squares must be used. One such method that is applicable to NARMAX model 
estimation is extended least-squares. This technique provides unbiased param- 
eter estimates and is iterative [19,21,22]. 

4.1 Extended Least-Squares 

Extended least-squares is one method, appropriate for NARMAX models that 
easily enables unbiased estimates to be computed. ELS addresses the bias prob- 
lem by modeling the lagged errors to obtain unbiased parameter estimates. 
ELS for linear systems has been widely studied and is also referred to as 
Panuska’s method, the extended matrix method, or approximate maximum 
likelihood [7-9]. 

In general, since the noise sequence is a realization of a stochastic process, 
it is not possible to solve for the noise source e, and it will not be equal to the 


prediction errors [23]. The prediction errors, e £ R iVx 1 , are defined as 

£ = Z - Z (5) 

where Z £ M iVx 1 is the measured output and Z = V P(9 £ M ,Vx 1 is the predicted 
output. In ELS, the NARMAX formulation of Eqn. 1 is redefined into a predic- 
tion error model by replacing e with £; making it a deterministic least-squares 
problem. 

The ELS formulation is defined as 

e = (^ T ^)~ l ^ T Z- where (6) 

9 £ M pxl is an ELS estimate of the system parameters, V P £ M ,Vx/ ' is a parti- 
tioned regressor matrix where l P, u is a function of z and u only, 1 P ZM g represents 
all the cross products involving £, and 'Pg is a polynomial function of the pre- 
diction errors only [23]. 

Often, ELS is considered a pseudolinear approach to parameter estimation 
[7,9,24]. Strictly speaking, the introduction of prediction errors into the model 
formulation no longer makes the model linear-in-the -parameters because the 
prediction errors depend on the model output, which is a function of all model 
parameters. The ELS technique solves a nonlinear optimization problem by 
ignoring the nonlinear character of the model and employing a least-squares 
approach. Essentially, ELS uses an approximate gradient of the model output 
with respect to the model parameters as a regression vector. 


5 Standardized Least-Squares 


In many nonlinear systems there are large numerical differences of the re- 
gressors due to the basis function(s) used to estimate the static map. The 
large differences lead to the regressor matrix being ill-conditioned and results 
in unstable matrix inversion and poor parameter estimates [25]. To alleviate 
ill-conditioning we propose using standardized least-squares (SLS). The SLS 
technique is summarized as follows. 

Given a matrix of independent variables V P and of dependent variables Z 
compute the mean and standard deviation of each variable and replace V P and 
Z with the centered and standardized variate as 


*P = OP - and Z = (Z-At z )Iz 1 (7) 


where £>i< is a diagonal matrix of standard deviations with denoting the 
standard deviation of the kth column of 'P and /./ ^ is a matrix whose kth column 
has all entries equal to the mean of column k of V P. Substituting Eqn. 7 into 6 
yields a SLS formulation. The SLS objective function is compactly expressed 
as 


1 

min - 
e 2 


(z-*pe) 


2 

2 


( 8 ) 


This is an unbiased estimator and converges asymptotically to the true system 
parameters. 
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6 Simulations 


The efficacy of the SLS ELS technique was assessed to provide stable estimates 
of continuous-time parameters describing ankle dynamics by simulating the 
parallel pathway model in continuous-time using Simulink (Fig. 1). The inputs 
were bandlimited 30 Hz pseudo-random binary sequences (PRBS) which were 

0.02 rad (peak-to-peak) and had a 35 ms switching rate. A PRBS input was 
used to excite the system dynamics since it is standard practice in experimental 
settings. The continuous-time coefficients used in this study correspond to 
values found in experiments and are shown in Table 2. One thousand Monte- 


CT Coefficient 

Value 

/ 

0.015 N m/rad/s 2 

B 

0.8 Nm/rad/s 

K 

150 Nm/rad 

CO 

40.0 

c 

1.00 

8 

10.00 Nm/rad/s 

A 

0.045 s 

T 

0.005 s 


Table 2. Continuous-time coefficient values. /: inertia, B: viscosity, K: elas- 
ticity, co: natural frequency, C: damping parameter, g: reflex stiffness gain, A: 
reflex delay and T : sampling interval. 

Carlo simulations were generated in which each input-output realization was 
unique, and had a unique Gaussian, white, zero-mean, noise sequence added to 
the output. Excluding the noise-free (NF) case the signal-to-noise ratio (SNR) 
of the noise sequence was decreased from 20 to 0 dB in increments of 5 dB. 
For identification, the data length was N = 4,000 points. 

6.1 Continuous- Time Parameter Estimation 

The system parameters were estimated as outlined in Eqns. 6-8. Specifically, 
for each input-output realization, we analyzed the ankle dynamics model as 
follows: 

1. Non- Standardized Extended Least-Squares : Ankle dynamics were ana- 
lyzed using the ELS estimator. 

2. Standardized Extended Least-Squares : Ankle dynamics were analyzed 
using the ELS estimator and standardized data. 

Fig. 2 shows a typical input-output realization used for this analysis. Employing 
our approach, continuous-time parameters were computed from discrete-time 
NARMAX estimates using the theoretical relationships given in Table 3. 
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PRBS Position Input 



Torque Output 



Figure 2. Typical PRBS ankle position input and torque output with 20 dB 
SNR used for Monte-Carlo analysis. 


CT Coefficient 

DT Relationship 

/ = 

B = 
K = 

CO = 

e = 

g = 

06+01 Y1 _ rll . B \ i 21 B T 

— V 7-2 T j ) T ji t 1 

TSgt" 1 - (o 2 7--) x J, 

- 2-202 _ 

— 1 + 02 — 0 i COT" 

08 x (4 + w 2 r 2 + 4£ft>r) x T = S^r 


Table 3. Discrete to continuous-time relationships for parameters I, B, K, 0 ), C 
and g of the parallel pathway ankle model. 
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The results of this study are presented in Fig. 3. The panels show plots of 
standard deviation about the mean for estimated continuous-time parameters of 
the linear and nonlinear paths (I, B, K, g, a> and C) using the non-standardized 
and standardized ELS identification techniques. Column one shows the linear 
path parameters. The continuous-time estimates using non-standardized ELS 
and relationships given in Table 3 provided a mean value that was close to 
the true parameter and had monotonically increasing variance from 20-10 dB 
SNR. Flowever, with decreasing SNR (5-0 dB) the parameter mean was biased 
away from the true with non-monotonically increasing variance that included 
estimates with incorrect sign. Using standardized data with the ELS algorithm, 
mean and variance estimates for the linear path parameters improved for all 
SNR levels but experienced a similar non-monotonically increasing variance 
for 5-0 dB SNR and encompassed values with opposite sign. Column two 
represents the nonlinear path estimates. Applying non-standardized and stan- 
dardized ELS to the nonlinear path illustrates that both approaches gave sim- 
ilar results. For the damping (£) and gain (g) parameters standardized data 
provided estimate that were worse than non-standardized data. The natural fre- 
quency (to) estimates were moderately improved using standardized data. Re- 
sults show that, for the linear path, the parameters were close to the true mean 
but variance estimate was non-monotonic and included estimates with incor- 
rect sign for low SNR levels. However, for the nonlinear path, continuous-time 
parameters computed using both approaches were significantly biased from the 
true mean for all SNR levels but did not include parameter estimates with in- 
correct sign. 


7 Discussion 

In this section, we summarize the findings of this study and discuss possible 
future approaches. The discussion highlights major features of our approach 
and offers possible explanations for the unexpected errors. The suggestion for 
future direction provided two broad alternatives to our solution, namely, an 
alternative identification strategy and discretization approach. 

7.1 Continuous- Time Parameter Estimation 

Estimation of continuous-time parameters of ankle dynamics demonstrates that, 
with noise-free input and output the both the ELS and SLS ELS approaches 
provided estimates that agree with the true known values. However, with out- 
put additive noise the linear path parameter estimates had non-monotonically 
increasing variance and encompassed values with opposite sign (see column 1 , 
Fig. 3). For the nonlinear path, the parameters were biased and the standard 
deviation did not encompass the true value (see column 2, Fig. 3). Although 
the nonlinear path parameters were biased, they did not experience sign flips as 
the linear path parameters did. 

The reason for the parameter bias may be due to the opposing nature of 
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Inertial Parameter 


Damping Parameter 






NF 20 15 10 5 


Signal-to-Noise Ratio (dB) 


„ „4Viscous Parameter 
x 10 



4 NF 20 15 10 5 0 

Signal-to-Noise Ratio (dB) 


. _6Elastic Stiffness 

x 10 


~ 5 

"D 

CO 



10 NF 20 15 10 5 0 

Signal-to-Noise Ratio (dB) 


Signal-to-Noise Ratio (dB) 

Natural Frequency 



Signal-to-Noise Ratio (dB) 

Gain 



Signal-to-Noise Ratio (dB) 


Figure 3. Non-standardized and standardized ELS approaches to continuous- 
time parameter estimation of ankle dynamics. Column 1: Linear path param- 
eters, 7, B,K. Column 2: Nonlinear path parameters, Ordinate: STD 

about mean. Abscissa: Output SNR = NL, 20, 15, 10, 5, 0 dB, where NL 
denotes noise-free case. (Note that the abscissa is shown in decreasing SNR 
which corresponds to increasing noise amplitude.) 
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the two paths. The linear path is described by a finite impulse response (FIR) 
system, which has theoretical poles close to the Nyquist frequency; requiring a 
high frequency signal to excite the dynamics sufficiently. The nonlinear path is 
described by a combination of a nonlinear static function and linear dynamics. 
The nonlinear function is described as a half-wave rectifier. This is considered a 
hard nonlinearity because it generates multiple harmonics and its output cannot 
be written analytically. To limit the effect of higher order harmonics and avoid 
internal aliasing at the nonlinearity output the input signal was bandlimited to 
be of low frequency. The opposing requirements of the two paths are central to 
this estimation problem and difficult to satisfy simultaneously. It is not possible 
to make more definitive comments about this approach without further detailed 
investigation and analysis. Below we outline a few alternative approaches to 
solve this identification problem. 

7.2 Closed-Loop Simulation 

A possible solution to continuous-time parameter estimation for ankle dynam- 
ics is to identify the model in closed-loop, removing the effects of derivatives 
(FIR realization of the linear path). In feedback, intrinsic stiffness is repre- 
sented as a compliance model (pure infinite impulse response (HR) realization). 
Although it may be possible to pose the identification problem as a feedback 
model, the true system is represented from position input to torque output. The 
recorded output torque has observation noise associated with it. When iden- 
tifying the model in feedback the noisy output is considered the input. Noise 
on the input violates assumptions and conditions for least-squares to yield un- 
biased parameter estimates. Implementing a feedback identification approach 
may require advanced methods such as total-least squares [26,27]. 

7.3 Discrete to Continuous- Time Parameter Mapping 

Another possibility for the source of error for continuous-time parameters com- 
puted using our NARMAX approach may be related to the nonlinear relation- 
ships between discrete-time NARMAX parameters and continuous-time pa- 
rameters of the physical system (see Table 3). A small deviation from the true 
parameter value in discrete-time, due to noise or numerical error, may appear 
as a significant error in the estimated continuous-time coefficients. As a result 
it may be advantageous to study these model parameters only in discrete-time 
as a combination of physiological effects. 

7.4 Continuous to Discrete-Time Transformations 

To derive a NARMAX model formulation the bilinear transform and Newton’s 
backwards formula were used to convert the continuous-time linear dynam- 
ics to discrete-time. These transforms were implemented since both require 
only a simple substitution to convert a continuous-time system to discrete-time. 
Two other techniques that give better approximation for linear time invariant 
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(LTI) systems are linear extrapolation and linear interpolation methods [28,29]. 
However, these techniques are seldom used due to added complexity for little 
gain. 

Linear extrapolation gives an improper transfer function (i.e. more zeros 
than poles) and linear interpolation gives a transfer function that is strictly 
proper (i.e. equal zeros and poles) [28,29]. These methods produce a discrete- 
time transfer function that gives a better output response than the bilinear trans- 
form or Newton’s backwards formula. However, the pitfall of these methods 
is that it is difficult to derive the coefficients of the discrete-time linear system. 
This is the main reason that almost all engineering text books and literature 
only discussed the bilinear transform and Newton’s backwards formula. 

Using linear interpolation or linear extrapolation to derive a NARMAX 
representation of ankle dynamics it may be possible to derive a better approx- 
imation to the derivative than the one used i.e. Newton’s backwards formula. 
This may give a better approximation to the derivative and provide simulation 
data that matches the continuous-time process better. However, one drawback 
is that it may give a discrete-time approximation to the derivative that is higher 
than first order thereby increasing the complexity of the NARMAX represen- 
tation [28,29]. Another limitation is that the continuous to discrete mapping 
of the continuous-time parameters will be more complicated since it involves 
exponential functions. This may result in more sensitivity for continuous-time 
parameter estimation since small deviations from the true value in discrete-time 
will result in large errors in continuous-time. 

8 Conclusions 

Clearly, much work remains to be done to resolve the possible source(s) of 
error for continuous-time parameter estimation for the model of ankle dynam- 
ics. Our results show that the proposed standardized extended least-squares 
method provided some improvement for estimating the continuous-time pa- 
rameters of the linear path. However, the parameter variance was too high 
because it included parameters with incorrect sign. For the nonlinear path, this 
technique did not offer any improvement over the non-standardized extended 
least-squares approach. A possible explanation for these results is due to op- 
posing constrains needed to properly excite the dynamics of both paths simul- 
taneously. The linear path requires a high frequency signal while the nonlinear 
path needs low frequency. It is unclear how to satisfy both concurrently. It 
is not possible to make more definitive remarks about the root of this error 
and remains an open research question. We have offered several suggestions 
for alternative ways forward to estimate the continuous-time parameters of the 
parallel pathway model of ankle dynamics. Of these, we deem that identifying 
this model in closed loop and implementing more advanced estimation tech- 
niques such as total least-squares is likely the best way forward. 
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